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Abstract;  Computer  simulation  of  three-dimensional  incompressible 
flow  is  of  interest  in  many  navigation,  coastal,  and  geophysical  applica¬ 
tions.  This  report  is  the  the  fifth  in  a  series  of  publications  that  documents 
research  and  development  on  a  state-of-the-art  computational  model¬ 
ing  capability  for  fully  three-dimensional  two-phase  fluid  flows  with  ves¬ 
sel/structure  interaction  in  complex  geometries  (Farthing  and  Kees,  2008; 
Kees  et  ah,  2008;  Farthing  and  Kees,  2009;  Kees  et  ah,  2009).  It  is  pri¬ 
marily  concerned  with  model  verification,  often  defined  as  “solving  the 
equations  right”  (Roache,  1998).  Model  verification  is  a  critical  step  on 
the  way  to  producing  reliable  numerical  models,  but  it  is  a  step  that  is  of¬ 
ten  neglected  (Oberkampf  and  Trucano,  2002).  Quantitative  and  qualita¬ 
tive  methods  for  verification  also  provide  metrics  for  evaluating  numerical 
methods  and  identifying  promising  lines  of  future  research. 

Fully-three  dimensional  flows  are  often  described  by  the  incompressible 
Navier-Stokes  (NS)  equations  or  related  model  equations  such  as  the 
Reynolds  Averaged  Navier  Stokes  (RANS)  equations  and  Two-Phase  Reynolds 
Averaged  Navier-Stokes  equations  (TPRANS).  We  will  describe  spatial  and 
temporal  discretization  methods  for  this  class  of  equations  and  test  prob¬ 
lems  for  evaluating  the  methods  and  implementations.  The  discretization 
methods  are  based  on  stabilized  continuous  Galerkin  methods  (variational 
multiscale  methods)  and  discontinous  Galerkin  methods.  The  test  prob¬ 
lems  are  taken  from  classical  fluid  mechanics  and  well-known  benchmarks 
for  incompressible  flow  codes  (Batchelor,  1967;  Chorin,  1968;  Schafer 
et  ah,  1996;  Williams  and  Baker,  1997;  John  et  ah,  2006).  We  demonstrate 
that  the  methods  described  herein  meet  three  minimal  requirements  for 
use  in  a  wide  variety  of  applications:  1)  they  apply  to  complex  geometries 
and  a  range  of  mesh  types;  2)  they  robustly  provide  accurate  results  over  a 
wide  range  of  flow  conditions;  and  3)  they  yield  qualitatively  correct  solu¬ 
tions,  in  particular  mass  and  volume  conserving  velocity  approximations. 


Disclaimer:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes.  Citation 
of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products.  All  product 
names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to  be  construed  as 
an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 

DESTROY  THIS  REPORT  WHEN  NO  LONGER  NEEDED.  DO  NOT  RETURN  IT  TO  THE  ORIGINATOR. 
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1  Introduction 


Fluid  flow  in  the  vicinity  of  vessels  and  structures  typically  becomes  quite 
complex  for  even  moderate  flow  conditions.  Flow  conditions  are  typically 
characterized  by  the  dimensionless  Reynolds  number  (Re)  given  by  Re 
=  VL/v  where  v  is  the  kinematic  viscosity  and  V  and  L  characterize  the 
velocity  and  length  scale  of  a  given  problem.  Figure  i  shows  a  von  Kar- 
man  vortex  street,  which  is  a  well-known  unsteady  flow  pattern  that  can 
develop  behind  a  cylinder  at  Re  near  too,  long  before  the  onset  of  tur¬ 
bulent  flow  (Batchelor,  1967).  For  open  channel  and  coastal  modeling, 
turbulent  flows  must  be  approximated  using  a  Reynolds  averaging  for¬ 
malism  (RANS)  or  Large  Eddy  Simulation  (LES).  Complex  averaged  ve¬ 
locities  develop  in  these  large  scale  flows  regardless  (Hutter  and  Johnk, 
2004).  Eluid  flows  with  a  free  surface,  such  as  two-phase  air /water  flow, 
introduce  significantly  more  complexity  because  the  motion  of  the  free 
surface  (waves)  induces  additional  velocity  variation  in  space  and  time. 
The  phenomenon  that  is  primarily  responsible  for  generating  temporal 
and  spatial  complexity  in  these  flows  is  the  interaction  of  the  strongly 
nonlinear  inertial  terms  with  the  weak,  small-scale  viscous  terms  in  the 
equations.  Numerical  methods  for  solving  this  class  of  equations  must 
address  the  destabilizing  influence  of  the  inertial  (advective)  terms  in  or¬ 
der  to  obtain  accurate  solutions. 

When  simulating  turbulent  flows,  two-phase  flows,  or  turbulent  two- 
phase  flows,  the  velocity  field  must  be  used  to  drive  additional  trans¬ 
port  equations  such  as  the  turbulence  closure  models  and  the  free  sur¬ 
face  models.  If  the  computed  velocity  field  does  not  satisfy  the  continuity 
equation,  then  this  error  leads  to  incorrect  results  in  these  models.  The 
error  then  propagates  to  the  turbulent  and  free  surface  models,  which  in 
turn  feeds  back  to  the  flow  model.  Attention  must  be  paid  to  the  “com¬ 
patibility”  of  numerical  solutions,  particularly  in  these  cases  of  coupled 
flow  and  transport  (Dawson  et  ah,  2004). 

In  this  report  we  describe  finite  element  discretizations  for  variable  coef¬ 
ficient  NS  equations  that  use  a  multiscale  approach  to  stabilizing  the  mo¬ 
mentum  advection  term.  The  approach  applies  to  unstructured  meshes 
and  variable  order  polynomial  approximation  spaces.  Eurthermore,  we 
employ  a  post-processing  approach  that  produces  locally  conservative  ve- 
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Figure  1.  Periodic  vortex  shedding  at  Re  100. 


locity  approximations  as  well  as  an  adaptive,  variable  order,  variable  step 
size,  temporal  discretization.  We  apply  the  methods  to  a  range  of  test 
problems  to  verify  the  implementation  and  evaluate  its  accuracy. 

The  outline  of  the  remainder  of  this  report  is  as  follows.  We  begin  by 
presenting  flow  formulations  representative  of  the  class  of  equations  for 
which  the  discretizations  are  applicable.  Then  we  present  details  on  the 
variational  multiscale  method  applied  to  a  representative  flow  model  as 
well  as  the  time  discretization  and  velocity  post-processing.  We  consider 
a  range  of  test  problems  to  verify  the  correctness  of  the  implementation 
and  analyze  the  results  to  achieve  a  better  understanding  of  the  robust¬ 
ness,  accuracy,  and  efficiency  of  the  mathematical  models.  We  conclude 
with  some  recommendations  for  future  research  and  development  on  nu¬ 
merical  methods  for  this  class  of  problems. 
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2  Formulation 


We  begin  with  a  physical  domain  O  and  a  time  interval  [0,  T].  We  write 
the  NS  equations  for  an  incompressible,  Newtonian  fluid  in  O  x  [0,  T]  as 

V  •  V  =  0  (1) 

9v  1 

—  +  V  •  Tv  0  V  -  v(Vv  +  VvOl  =  g - Vp  (2) 

ot  ^  ^  p 

where  v  is  the  velocity,  v  0  v  is  the  tensor  {niUj} ,  i,j  =  1, 2, 3,  v  is  the 
kinematic  viscosity,  g  is  the  gravitational  acceleration,  and  p  is  the  den¬ 
sity.  Fully  describing  either  RANS  or  TPRANS  models  is  beyond  the 
scope  of  this  report  and  we  merely  give  a  representative  formulation.  A 
RANS  formulation  with  a  first  order  turbulence  closure  model  can  be 
written  as 


V  •  V  =  0  (3) 

+  V  •  [v0v- (v  + Vf)(Vv  + V^)]  =  g--V(p  +  ^)  (4) 

ot  ^  ^  p  3 

where  v  and  p  are  the  Reynold’s  averaged  velocity  and  pressure,  Vf  is  the 
turbulent  kinematic  viscosity,  and  k  is  the  turbulent  kinetic  energy  (Hut- 
ter  and  Johnk,  2004;  Bernard  et  ah,  2007).  In  this  case  v  and  p  are  the 
unknowns  and  Vt  and  k  are  also  part  of  the  solution  arising  through  the 
coupling  to  a  turbulence  closure  model.  For  an  air/water  flow,  neglecting 
the  effect  of  surface  tension,  we  can  write  the  TPRANS  model  equations 
as 


V  •  V  =  0  (5) 

^  +  V  •  [v  0  V  -  (vfcp)  +  vd(p))(Vv  +  V^)]  =  g  -  +  y ) 

where  (/>  is  a  function  describing  the  fluid  distribution  (e.g.  a  level  set  or 
volume  of  fluid  function).  In  this  case  q)  is  an  additional  solution  vari¬ 
able  arising  through  the  coupling  of  an  equation  for  the  fluid-fluid  inter¬ 
face. 
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2.1  Test  equation 

Since  our  focus  in  this  report  is  specifically  on  spatial  discretizations  for 
the  flow  equation  and  not  on  turbulence  or  free  surface  modeling,  we 
will  focus  on  the  general  variable  coefficient  NS  equation 

V  •  V  =  0  (7) 

9v  1 

—  +  V  •  [v  0  V  -  v(x)(Vv  +  VvO]  =  g - 

dt  ^  ^  p(x) 

Henceforth  we  will  drop  the  explicit  dependence  on  x. 

2.2  Weak  formulation 


We  proceed  by  defining  a  standard  weak  formulation  of  the  NS  equation. 
Boundary  conditions  are  an  important  and  complex  aspect  of  real  world 
modeling  that  we  will  not  treat  fully  in  this  report.  Instead  we  will  as¬ 
sume  that  the  boundary  of  the  domain  has  two  partitionings:  [dY^,  9r^} 
and  Y'^}  on  which  the  boundary  conditions  are  given  as 


p  =  Pd  on  Ff, 

(9) 

V  •  n  =  on  F^ 

(10) 

V  =  Vp  on  F^ 

(11) 

v(x)(Vv  +  VvO]  •  n  =  on  F](r 

(12) 

Furthermore  we  assume  that  p(x,  0)  =  po  and  v(x,  0)  =  Vq  are  pre¬ 
scribed  initial  conditions.  Since  our  focus  is  on  numerical  methods  and 
test  problems,  we  state  an  abstract  weak  formulation  of  NS  problems 
leaving  out  almost  all  rigorous  details  except  those  necessary  to  define 
the  numerical  methods.  First,  we  will  seek  a  solutions  p  and  v  that  are 
members  of  spaces  of  functions  Vf(0,  T;  W(0)),  and  Vj(0,  T;  V'^(O)).  In 
particular  this  means  thatp(f)  G  W(0),  v(f)  g  V^(0)  and  that  the  Dirich- 
let  boundary  conditions  are  incorporated  into  the  definition  of  W(0) 
and  V^(0).  We  say  a  solution  is  a  weak  solution  if 


-  /  V  •  Vw^dV 

Ja 


9v 


/  -  V  0  V  •  Vw^dV 

la  dt 


hPdS  \/wP  G  W^(0) 


v(x)(Vv  +  W)  ■  Vw'^dV 


(13) 


Ja 

^  la  P(x) 


Vp  w'^dV 


Vw"  G  W^(0) 


lyu 


(14) 
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where  we  interpret  vector-vector  multiplication  as  component-wise  mul¬ 
tiplication  (i.e.  equation  14  is  a  vector  equation).  We  call  W(0)  the  trial 
space  forp  and  the  test  space  forp. 
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3  Discrete  Approximation 


We  now  define  finite  dimensional  approximation  spaces  corresponding 
to  the  abstract  function  spaces  above.  This  converts  that  abstract  weak 
formulation  into  a  problem  on  the  finite  dimensional  vector  space 
where  N  is  the  number  of  discrete  degrees  of  freedom. 

3.1  Time  discretization 


First  we  partition  the  time  interval  as  [to,  tn,  tn+i,  •  •  • ,  T].  Our  choice 

of  space  for  V^fO,  T;  V^(0))  will  be  a  subspace  of  the  continuous  func¬ 
tions,  Co(0,  T ;  V^(0)),  including  certain  polynomials  defined  on  the  time 
discretization.  In  particular,  we  will  assume  that  for  n  +  1  >  /c  that  u  is  a 
Lagrange  polynomial  in  t  of  the  form 

rik 

v(f)  =  ^  x)  (15) 

k=0 

where  4  is  the  Lagrange  basis  function  at  tn+i-k-  Assuming  vitn+i-k)  is 
known  for  k  >  0,  this  implies  that 


— (4+i)  =  av(f„+i,x)  +  p 
ot 

(16) 

where  a  and  P  depends  on  {4}  and  {^itn+i-k,  ^)}  for  k  =  1, . 
To  simplify  the  notation  we  define 

..,nk  and  4. 

DfVn+i(x)  :=  av(f„+i,x)  +  p 

(17) 

This  approximation  converts  the  initial-boundary  value  problem  into  a 
sequence  of  boundary  value  problems  at  ti,  Dropping  the  time 

subscript  n  +  1  we  write  the  weak  formulation  of  the  boundary  value 
problem  as 


v-S/ufdV  = 


la. 


hPdS  Wuf  e 


(18) 


DfVW  -  V  (8)  V  •  Vw^dV  = 


la 


v(x)(Vv  +  W)  ■  Vw^dV 


Ja 

^  la  P(x) 


Vp  w'^dV 


h^^dS  Vw^  e 


(19) 
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3.2  Multiscale  formulation 

We  now  build  an  approximate  weak  formulation  in  time  using  the  mul¬ 
tiscale  formalism  of  (Hughes,  1995).  Let  M*  be  a  simplicial  mesh  on  O 
in  Ud  =  2, 3,  containing  Ng  elements,  {Og},  e  =  1, .  ..,Ne,  Nf  faces, 
{y/},/  =  1, .  ..,Nf,  and  Nn  nodes,  {x^},  n  =  1, .  ..,Nn.  The  collection  of 
faces  in  the  domain  interior  is  denoted  F/.  We  also  assume  that  the  in¬ 
tersection  of  elements  Oe,  Og'  G  is  either  empty,  a  unique  y/  G  F/, 
an  edge  (for  R^),  or  a  point.  The  diameter  of  Og  is  hg  and  its  unit  outer 
normal  is  written  Ug. 

Consider  test  and  trial  spaces  V  and  W.  The  basic  idea  of  a  multiscale 
method  is  to  split  V  and  W  into  resolved  and  unresolved  scales  using 
direct  sum  decompositions 


=  Vh®8V 
W  =  Wh®6W 


(20) 


(21) 


For  this  work  Vh  and  Wh  are  the  continuous,  piecewise  polynomial  spaces 
of  the  classical  Galerkin  finite  element  method: 


Vh  =  {vheVnC%ay.Vh\n,eP\ng)} 
Wh  =  {lOhG  WnC°(0):mh|a,  GP'^(Og)} 


(22) 


(23) 


while  8V  and  8W  remain  infinite  dimensional.  We  will  consider  k  =  1 
or  k  =  2  and  the  standard  extensions  of  these  spaces  to  spaces  of  vector 
valued  functions  written  as  V/,  and  W^. 

With  this  decomposition  for  W  and  W,  the  solution  is  written  uniquely 
as  p  =  Ph+p'  and  v  =  v/,  +  v'.  After  some  manipulation  and  approximation 
(Hughes,  1995),  we  obtain  the  weak  formulation 


(25) 


L 


AT 


vwj;  G  wi;(n) 
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where 


(26) 


(28) 


(27) 


The  operators  and  L*p^^  are  the  adjoint  operators  corresponding  to 
the  divergence  and  pressure  gradient  operators.  The  operator  L*  is  the 
adjoint  of  the  operator  obtained  by  linearizing  the  first  term  in  25  and 
assuming  that  DtWh  and  V  •  v  are  zero. 

3.3  Algebraic  sub-grid  scale  approximation 

To  obtain  a  closed  set  of  equations  for  ph,  we  need  approximations  for 
p'  and  v'.  We  use  the  standard  Algebraic-Sub-Grid  Scale  (ASGS)  approxi¬ 
mations  given  by 


p'  =  -TpRp 


(29) 


(30) 


where 


4v  +  2p||Vft,„_i||/i+  \D'^\h^ 


(31) 


1 


^  ,  2p||Vft,n-i|| 

h2  h 


(32) 


and 


(33) 


(34) 


(35) 


3.4  Velocity  post-processing 


Quite  often  a  velocity  approximation  along  the  boundaries  of  the  mesh 
elements  is  required  as  input  to  other  models  such  as  chemical  species 
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transport  and  particle  tracking.  One  shortcoming  of  the  finite  element 
approximations  above  is  that  the  condition  V  •  v  =  0  may  not  be  locally 
conservative,  i.e. 


V/,  •  ndS  ^  0 


(36) 


For  this  reason  we  post-process  v  to  obtain  a  new  velocity  v  in  the  space 
V  defined  by 


V(0)  =  {v  G  Co(0)  :  via,  G  ©  (xP°(0,))}  (37) 

The  space  V(0)  is  the  velocity  space  for  the  well-known  Raviart-Thomas 
space  of  order  zero  (RTo).  The  post-processed  velocity  satisfies  equation 
36  up  to  the  accuracy  of  the  nonlinear  solver  and  the  error  is 

||v- V*||(L2(fl))"d  <  C/l  (38) 

for  some  constant  C  depending  on  the  exact  solution  v*  but  independent 
of  the  maximum  element  diameter  h.  This  approach  was  originally  pre¬ 
sented  in  (Larson  and  Niklasson,  2004)  and  the  implementation  in  this 
work  was  evaluated  on  a  variety  of  variable  coefficient  test  problems  and 
unstructured  meshes  in  (Kees  et  ah,  2008). 

3.5  Additional  details 

The  discretization  above  yields  a  system  of  nonlinear  algebraic  equations 
at  each  time  step.  To  solve  these  systems  in  the  time-dependent  case  we 
used  Newton’s  method.  In  steady  state  cases  we  used  either  Newton’s 
method  or  Pseudo-transience  continuation  (Knoll  and  McHugh,  1998; 
Farthing  et  ah,  2003).  The  2D  simulations  were  run  on  a  MacPro  2x3 
GHz  Quad-Core  Intel  Xeon  processor  with  16  GB  of  memory.  On  this 
system  linear  systems  were  solved  using  the  SuperLU  (serial)  sparse  di¬ 
rect  solver  (Demmel  et  ah,  1999).  The  3D  simulations  were  run  on  32 
or  64  cores  of  a  Dell  Linux  Cluster  with  1955  2x2.66  GHz  Quad-Core 
Intel  Xeon  nodes  with  8GB  of  memory.  On  this  platform  we  used  the 
SPOOLES  parallel  sparse  direct  solver  (Ashcraft  and  Grimes,  1999)  via 
the  PETSc  framework  (Balay  et  ah,  2001,  2004,  1997). 
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4  Model  Verification 

4.1  Plane  Poiseuille  and  Couette  Flow 

First  we  consider  steady-state  flow  between  two  parallel  plates  of  infinite 
extent,  where  the  flow  is  driven  by  the  movement  of  the  top  plate  and/or 
an  externally  applied  pressure  gradient.  Assuming  the  Z-axis  is  normal 
to  the  plates  and  that  the  flow  and  pressure  gradient  are  alligned  with 
the  Z-axis,  the  solution  to  the  incompressible  NS  equations  in  this  case 
is 


7TT  1  Fin* 

u\Z)  =  +  ^ZiZ-Lz)  V=W  =  0  (39) 

Lz  2fi  dX 

P*(X)  =  ^X+po  (40) 

where  H  is  the  distance  between  the  plates,  (U,  0, 0)  is  the  velocity  of  the 
top  plate  relative  to  the  bottom  plate  and  po  is  an  arbitrary  constant.  To 
use  this  solution  for  model  verification  on  a  finite  domain  we  consider 
a  translated  and  rotated  coordinate  system  (x,  y,  z)  and  a  rectangular  re¬ 
gion  between  the  plates  given  by  O  =  [0, Lx]^[0, Ly]x[0, L^].  In  particular 
we  use 


X  =  X  •  n„  (41) 

Z  =  Xtlp-Zs  (42) 

tip  =  (cos(0p)  sin(^p),  sin(0p)  sin(^p),  cos(^p))^  (43) 

tiu  =  (cos(0„)  sin((pi,),  sin(0u)  sin((/)u),  cos((pu))^  (44) 


To  verify  the  spatial  discretizations  we  solved  this  problem  with  four  lev¬ 
els  of  mesh  refinement  choosing  Op,  cpp,  O^,  (pv  so  that  flow  is  skew  to  the 
grid.  The  results  for  2D  with 


<Pp 

e. 

(pv 

Zs 

2jr/3 

jt/2 

Ji/Q 

jt/2 

-1/2 

and  3D  with 
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Table  1.  Grid  refinement  study  for  2D  Poiseuille  problem. 


h 

||p-p1|L2(n) 

||u-u*||L,(n) 

\\v-v*\\L,in) 

dX  ■ 

II 

o' 

II 

0.5 

1.53  X  10“® 

5.82x10-12 

7.52x10-12 

0.25 

9.66x10-“ 

7.03x10-1® 

1.03x10-12 

0.125 

6.97x10-11 

1.12x10-1® 

4.80x10-14 

0.0625 

5.17x10-12 

8.43x10-1® 

7.65x10-1® 

II 

-1,A:=  1 

0.5 

1.01  X  10^ 

1.17  X  lOi 

6.71  X  10° 

0.25 

1.16  X  10^ 

3.09  X  10° 

1.69  X  10° 

0.125 

3.25  X  lO"! 

7.37  X  10-1 

4.64  X  10-1 

0.0625 

8.97  X  lO^ 

1.78  X  10-1 

1.29  X  10-1 

dp  . 
dX  ■ 

CM 

II 

o' 

II 

0.5 

2.77  X  10-® 

2.33x10-12 

2.64x10-12 

0.25 

4.02x10-1° 

7.81x10-1® 

4.29x10-1® 

0.125 

9.58x10-11 

1.21x10-1® 

8.85x10-14 

0.0625 

1.16x10-11 

2.03x10-1'^ 

1.10x10-14 

II 

I 

II 

0.5 

6.66  X  10-° 

4.10x10-14 

2.88x10-14 

0.25 

5.49  X  10-° 

5.34x10-14 

2.61x10-14 

0.125 

1.11x10-® 

8.80x10-14 

6.74x10-14 

0.0625 

2.12x10-® 

1.72x10-1® 

1.33x10-1® 

Op 

<Pp 

0, 

(pv 

Zs 

Ji/Q 

Ji/Q 

Ji/3 

Ji/6 

0 

are  given  in  Figures  i  and  2.  The  quadratic  finite  element  approximation 
(k  =  2)  is  accurate  to  within  the  nonlinear  solver  tolerance  of  1.0  x  10“^ 
in  all  cases,  demonstrating  that  it  is  essentially  able  to  represent  the  true 
solution  exactly.  The  linear  finite  element  approximation  (k  =  1)  is  es¬ 
sentially  exact  when  the  solution  is  linear  (^  =  0)  and  demonstrates 
quadratic  convergence,  which  is  consistent  with  the  theoretical  a  priori 
error  estimates  for  smooth  solutions. 

4.2  Vortex  decay 

This  time-dependent  problem  was  originally  described  in  (Chorin,  1968) 
and  was  used  to  study  time  discretizations  for  the  NS  equations  in  (John 
et  ah,  2006).  The  flow  domain  is  again  O  =  [0, 1]  x  [0, 1].  The  analytical 
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Table  2.  Grid  refinement  study  for  3D  Poiseuille  problem. 


h 

\\p-p*\\L,m 

l|u-i^1U2(n) 

l|u-u1U2(n) 

\\w-w*\\L,(n) 

^  -  0  k  - 

1 

0.5 

3.28  X  10“^ 

3.94x10““ 

1.92x10““ 

6.30x10““ 

0.25 

4.34x10““ 

1.18x10““ 

1.90x10““ 

4.21x10““ 

0.125 

1.18x10““ 

2.04x10““ 

2.61x10““ 

1.89x10““ 

0.0625 

3.12x10“^^ 

9.59x10““ 

1.47x10““ 

7.21x10““ 

1 

0.5 

4.02  X  10^ 

1.69  X  10^ 

2.91  X  10! 

1.75  X  10! 

0.25 

2.93  X  10^ 

4.55  X  10° 

7.70  X  10° 

4.87  X  10° 

0.125 

9.51  X  10"^ 

1.12  X  10° 

1.90  X  10° 

1.41  X  10° 

0.0625 

2.66  X  10“^ 

2.80x10“! 

4.70x10“! 

3.85  X  10“! 

2 

0.5 

8.89  X  10“^ 

1.95x10““ 

3.87x10“!! 

3.87x10“!! 

0.25 

2.74  X  10“3 

8.66x10““ 

1.53x10“!! 

1.10x10“!! 

0.125 

7.69x10““ 

1.32x10““ 

3.18x10““ 

2.10x10““ 

0.0625 

1.08x10““ 

3.05x10““ 

4.98x10““ 

3.34x10““ 

2 

0.5 

5.80  X  10“3 

3.41x10““ 

5.90x10““ 

4.79x10““ 

0.25 

9.92  X  10“^ 

4.40x10““ 

7.18x10““ 

6.30x10““ 

0.125 

2.74  X  10“® 

1.06x10““ 

2.80x10““ 

1.12x10““ 

0.0625 

2.11x10“® 

1.12x10““ 

2.19x10““ 

1.11x10““ 

solution  is  given  by 


P* 

=  --  (cos(2nu7rx)  sin(2nujry))  exp(-4n^jr^t/Re) 

(45) 

u* 

=  -  cos(ni,jTx)  sin(n:jjTy)  exp(-2n^jT^t/Re) 

(46) 

V* 

=  sin(ni,jTx)  cos(n„7ry)  exp(-2n^7r^f/Re) 

(47) 

and  Re  =  1/v.  We  use  this  solution  to  provide  non-homogeneous  Dirich- 
let  boundary  conditions  and  initial  conditions  for  all  variables.  The  so¬ 
lution  is  an  array  of  Uy  x  Uy  vortices  with  alternating  rotation  which  de¬ 
cay  in  time  exponentially  at  a  rate  controlled  by  Re.  In  Tables  3  and  4 
we  present  errors  for  refinement  in  space  and  time  for2  x  2  vortices  at 
Re  =  1  and  1  x  10®.  In  a  three  cases  the  high  Re  runs  failed  (denoted 
by  an  X)  due  to  repeated  reduction  of  the  time  step  in  cases  with  tol  < 
1.0  X  10““^.  This  failure  mode  occurs  because  the  ASGS  approximation 
becomes  badly  scaled  for  small  time  steps  (Bazilevs  et  ah,  2007).  Sub¬ 
grid  error  approximations  that  are  valid  for  small  timesteps  is  an  open 
area  of  research,  but  safeguarding  against  small  time  steps  (or  choos¬ 
ing  temporal  error  tolerances  appropriate  for  the  given  mesh)  should  be 
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sufficient  for  most  applications  where  a  time  step  on  the  order  of  the  ad- 
vective  Courant-Friedrich-Levy  condition  is  appropriate. 

4.3  Lid  driven  cavity 

The  flow  domain  is  O  =  [0,  a]  x  [0, 5]  x  [0,  c].  The  boundary  conditions 
are  given  by 

V  =  {U,  V", 0)  on  z  =  c,0  <  X  <  a,0  <  y  <  b 

V  =  0  on  X  =  0,a,y  =  0,b,z  =  0  (48) 

p{a/2,b/2,c)  =  0 

This  problem  has  no  analytical  solution  and  exhibits  a  wide  range  of  be¬ 
havior  depending  on  the  Re.  There  is  a  discontinuity  in  the  velocity  at 
the  boundary  along  the  upper  edges  of  the  cavity  (corners  in  2D).  The 
discontinuity  reduces  the  regularity  of  the  solution  and  consequently 
produces  a  reduction  in  the  asymptotic  order  of  convergence.  Never¬ 
theless,  it  is  a  standard  verification  problem,  and  a  great  deal  is  known 
about  the  structure  of  solutions  (Bassi  et  ah,  2006;  Erturk  et  ah,  2005). 
In  Figure  5  we  present  the  results  of  a  grid  refinement  study  with  four 
levels  of  mesh  refinement  using  the  fourth  level  as  the  “exact”  solution. 
This  measure  of  error  is  not  accurate  enough  to  compare  two  methods 
of  different  orders.  The  L2  error  estimates  are  shown  to  be  decreasing 
monotonically  but  clearly  the  order  of  convergence  is  less  than  the 
quadratic  and  cubic  rates  predicted  by  the  theory  for  smooth  solutions. 
The  streamlines  for  the  driven  cavity  in  2D  and  3D  are  given  in  Fig¬ 
ures  2-4.  The  structure  of  the  flows  is  in  close  agreement  with  previous 
numerical  studies,  in  particular  the  detailed  high  Re  studies  in  (Erturk 
et  ah,  2005). 

4.4  Backward  facing  step 

In  this  problem  we  consider  flow  over  a  square  step  at  the  lower  left 
hand  edge  of  the  domain.  We  can  describe  the  step  as 


Os  =  {x  :  0  <  X  <  as,  0  <  y  <  6, 0  <  z  <  Cs} 


(49) 
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Figure  2.  Lid  driven  cavity  in  2D  at  Re  =  100. 


the  flow  domain  is  then  O  =  [0,  a]  x  [0,  b]  x  [0,  c]  The  boundary 
conditions  are  given  by 


u 

V 

V 

V  •  n 

v(x,0,z) 
pix,  0,  z) 
P 

-v(x)(Vv  +  VvO  •  n 


z(z  -  on  X  =  0 

IV  =  0  on  X  =  0 


0 

0 

v(x,6,z) 

pix,b,z) 

0 

0 


on  90s,  z  =  0,  c 


on  90s,  z  =  0,  c 


on  X  =  a 
on  X  =  a 


(50) 


This  problem  has  been  the  focus  of  much  experimental  and  numerical 
study,  and  we  will  derive  the  details  of  the  test  problem  from  the  exper¬ 
imental  work  in  (Armaly  et  al.,  1983)  and  one  of  the  subsequent  numer¬ 
ical  studies  (Williams  and  Baker,  1997).  The  dimensions  of  the  domain 
are  given  by 


as 

Cs 

a 

c 

5 

4-9 

155 

10.1 
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Figure  3.  Lid  driven  cavity  in  2D  at  Re  =  20000. 


The  recirculation  length  of  the  primary  vortex  for  Re=  100  -  800  is 
given  in  Figure  5,  which  is  in  close  agreement  with  (Williams  and  Baker, 
1997).  Examples  of  the  vortex  structure  in  2D  and  3D  are  given  in  Fig¬ 
ures  6  and  7. 

4.5  Flow  past  a  cylinder 

We  consider  flow  around  a  cylinder  of  radius  R,  oriented  along  the  y- 
axis.  The  cylinder  can  be  described  implicitly  by 

Os  =  |x  :  <  r|  (51) 

the  flow  domain  is  the  O  =  [0,  a]  x  [0,  b]  x  [0,  c]  P|  where  is  the 
complement  of  Os.  The  boundary  conditions  are  given  by 


u 

on 

X  = 

0 

u 

=  v  =  0 

on 

X  = 

0 

V 

=  0 

on 

90, 

= 

0,  c 

V 

•0 

=  0 

on 

y  = 

0, 

b 

v(x)(Vv  +  VvO 

•  n 

=  0 

on 

X  = 

a, 

y 

=  0,b 

P 

=  0 

on 

X  = 

a 

This  problem  has  no  analytical  solution  and  exhibits  a  wide  range  of  be¬ 
havior  depending  on  Re.  The  variation  in  the  lift  coefficient  is  shown  in 
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Figure  8,  which  closely  matches  prior  results  studying  higher-order  time 
discretizations  (John  et  ah,  2006). 
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Table  3.  Grid  refinement  study  for  vortex  decy  probiem  Re  =  1. 


tol 

1.0  X  10-2 

1.0  X  10-2 

1.0  X  10-"^ 

1.0  X  10-2 

h  =  0.1, k  =  1 

\\p-p*\\L2m 

2.01x10-^ 

2.49  X  10-1 

2.59  X  10-^ 

2.54x10-^ 

\\u-u*\\L2m 

4.59  X  10-2 

5.64  X  10-2 

5.96  X  10-2 

5.94  X  10-2 

-  i^1lL2(n) 

4.59  X  10-2 

5.63  X  10-2 

5.95  X  10-2 

5.94  X  10-2 

/] 

=  0.05,  k=  1 

\\p-p*\\L2(n) 

1.22  X  10-^ 

8.29  X  10-2 

1.83  X  10-^ 

1.15x10-^ 

\\u-u*\\L2in) 

1.09  X  10-2 

1.19  X  10-2 

1.74  X  10-2 

1.72  X  10-2 

\\v-v*\\L2(n) 

1.09  X  10-2 

1.19  X  10-2 

1.73  X  10-2 

1.72  X  10-2 

h 

=  0.025,  k=  3 

llp-plli^cn) 

2.33  X  10-^ 

5.28x10-2 

4.45  X  10-2 

3.84x10-2 

\\u-u*\\L2(n) 

1.10x10-2 

3.14  X  10-2 

4.37  X  10-2 

4.39  X  10-2 

-  i’1lL2(n) 

1.10x10-2 

3.14  X  10-2 

4.37  X  10-2 

4.39  X  10-2 

h 

=  0.0125,  A:  = 

1 

||p-p1|L2(n) 

3.09  X  10-1 

2.79  X  10-2 

1.74  X  10-2 

1.19x10-2 

\\u-u*\\L2m 

1.29  X  10-2 

2.46  X  10-2 

1.03  X  10-2 

1.06  X  10-2 

-  i^1lL2(n) 

1.29  X  10-2 

2.46  X  10-2 

1.03  X  10-2 

1.06  X  10-2 

h  =  0.1, k  =  2 

\\p-p*\\L2m 

1.49  X  10-1 

8.68  X  10-2 

6.90  X  10-2 

5.85  X  10-2 

l|u-i^1lL2(n) 

1.33  X  10-2 

7.48  X  10-2 

1.22  X  10-2 

1.27  X  10-2 

1.33  X  10-2 

7.49  X  10-2 

1.23  X  10-2 

1.27  X  10-2 

/] 

=  0.05,  k  =  2 

llp-plU^Cn) 

2.41  X  10-1 

1.88  X  10-2 

1.52  X  10-2 

1.70x10-2 

llu-i^1lL2(n) 

1.36  X  10-2 

5.48  X  10-2 

1.00  X  10-2 

1.59  X  lO-'^ 

-  i^1lL2(n) 

1.36  X  10-2 

5.48  X  10-2 

1.00  X  10-2 

1.58  X  10-"^ 

h 

=  0.025,  A:  =  S 

\\p-p*\\L2(n) 

3.07  X  10-1 

1.88  X  10-2 

7.23  X  10-2 

3.66  X  10-2 

liu-  u*||L2(n) 

1.36  X  10-2 

6.53  X  10-2 

8.94  X  10-2 

3.84x10-2 

-  i^1lL2(n) 

1.36  X  10-2 

6.53  X  10-2 

8.92  X  10-2 

3.80  X  10-2 

h 

=  0.0125,  A:  = 

2 

||p-p1|L2(n) 

3.43  X  10-1 

1.51  X  10-2 

9.08  X  10-2 

9.15x10-'^ 

\\u-u*\\L2(n) 

1.36  X  10-2 

6.85  X  10-2 

9.52x10-2 

4.13x10-2 

||i^  - 

1.36  X  10-2 

6.85  X  10-2 

9.51  X  10-2 

4.13x10-2 
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Table  4.  Grid  refinement  study  for  vortex  decy  probiem  Re  =  1x10^. 


tol 

1.0  X  10-2 

1.0  X  10-2 

1.0  X  10-"^ 

1.0  X  10-® 

h  =  0.1, k  =  1 

\\p-p*\\L2m 

1.74  X  10^ 

2.24  X  10® 

2.33  X  10® 

2.30  X  10® 

\\u-u*\\L2m 

4.40  X  10-2 

5.44  X  10-2 

5.76  X  10-2 

5.74  X  10-2 

-  i^1lL2(n) 

4.40  X  10-2 

5.44  X  10-2 

5.76  X  10-2 

5.74  X  10-2 

/] 

=  0.05,  k=  1 

||p-p1|L2(n) 

1.16  X  10^ 

7.54  X  10"^ 

X 

1.09  X  10® 

\\u-u*\\L2in) 

1.07  X  10-2 

1.15  X  10-2 

X 

1.68  X  10-2 

\\v-v*\\L2(n) 

1.07  X  10-2 

1.15  X  10-2 

X 

1.68  X  10-2 

h 

=  0.025,  k=  3 

llp-plli^cn) 

2.33  X  10^ 

5.15  X  10^ 

4.35  X  10^ 

3.68  X  104 

\\u-u*\\L2(n) 

1.10x10-2 

3.10  X  10-® 

4.31  X  10-® 

4.33  X  10-® 

-  i’1lL2(n) 

1.10x10-2 

3.10  X  10-® 

4.31  X  10-® 

4.33  X  10-® 

h 

=  0.0125,  A:  = 

1 

||p-p1|L2(n) 

3.08  X  10^ 

2.76  X  10^ 

1.74  X  10^ 

1.01  X  104 

\\u-u*\\L2m 

1.29  X  10-2 

2.47  X  10-® 

1.03  X  10-® 

1.05  X  10-® 

-  i^1lL2(n) 

1.29  X  10-2 

2.47  X  10-® 

1.03  X  10-® 

1.05  X  10-® 

h  =  0.1, k  =  2 

llp-plU^Cn) 

1.51  X  10^ 

8.70  X  10^ 

5.77  X  10^ 

X 

|ju-  u*||L2(n) 

1.34  X  10-2 

7.51  X  10-® 

1.22  X  10-® 

X 

1.34  X  10-2 

7.51x10-® 

1.22  X  10-® 

X 

/] 

=  0.05,  k  =  2 

llp-plU^Cn) 

2.42  X  10® 

1.69  X  lO^ 

1.46  X  10^ 

X 

llu-i^1lL2(n) 

1.36  X  10-2 

5.31  X  10-® 

9.70  X  10-4 

X 

-  i^1lL2(n) 

1.36  X  10-2 

5.31  X  10-® 

9.70  X  10-4 

X 

h 

=  0.025,  A:  =  S 

||p-p1|L2(n) 

3.07  X  10® 

1.80  X  10"^ 

7.06  X  10® 

3.50  X  10® 

liu-  u*||L2(n) 

1.36  X  10-2 

6.51  X  10-® 

8.83  X  10-® 

4.12x10-® 

-  i^1lL2(n) 

1.36  X  10-2 

6.51  X  10-® 

8.83  X  10-® 

4.12x10-® 

h 

=  0.0125,  A:  = 

2 

||p-p1|L2(n) 

3.43  X  10® 

1.45  X  10"^ 

9.20  X  10® 

8.65  X  102 

\\u-u*\\L2(n) 

1.36  X  10-2 

6.86  X  10-® 

9.54x10-® 

4.13x10-® 

-  i’1lL2(n) 

1.36  X  10-2 

6.86  X  10-® 

9.54  X 10-® 

4.13x10-® 
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Table  5.  2D  lid  driven  cavity  at  Re  =  100,400,1000. 


h 

Wp-p'^h^m 

llu-u'^llL^cn) 

Re= 

100,A:  =  1 

0.1 

3.41  X  10-3 

5.85  X  10-3 

5.53  X  10-3 

0.05 

2.52  X  10-3 

3.32  X  10-3 

3.88  X  10-3 

0.025 

1.52  X  10-3 

1.39  X  10-3 

1.66  X  10-3 

Re= 

100,fc  =  2 

0.1 

1.16x10-3 

2.97  X  10-3 

3.26  X  10-3 

0.05 

7.96  X  10-“^ 

1.39  X  10-3 

1.56  X  10-3 

0.025 

4.98  X  10-“^ 

6.16  X  10"^ 

7.74  X  lO-'^ 

Re= 

400,A:  =  1 

0.1 

1.44  X  10-3 

2.88  X  10-3 

2.80  X  10-3 

0.05 

1.01  X  10-3 

1.52  X  10-3 

1.64x10-3 

0.025 

6.07  X  10-3 

5.27  X  10-3 

6.46  X  10-3 

Re= 

400,A:  =  2 

0.1 

6.71  X  10-3 

2.07  X  10-3 

1.82  X  10-3 

0.05 

3.88  X  10-3 

9.40  X  10-3 

8.73  X  10-3 

0.025 

2.12x10-3 

3.04  X  10-3 

3.40  X  10-3 

Re=  1000, A:  =  1 

0.1 

4.59  X  10-3 

9.04  X  10-3 

8.44  X  10-3 

0.05 

3.61  X  10-3 

6.29  X  10-3 

5.97  X  10-3 

0.025 

1.94  X  10-3 

2.53  X  10-3 

2.55  X  10-3 

Re=  1000, A:  =  2 

0.1 

4.34  X  10-3 

8.99  X  10-3 

8.22  X  10-3 

0.05 

2.87  X  10-3 

5.89  X  10-3 

5.32  X  10-3 

0.025 

1.18x10-3 

2.15x10-3 

1.99  X  10-3 

Figure  6.  Backward  facing  step  in  2D  at  Re  =  400  and  800. 
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velocity 

3.3e-04 

2.5e-04 

1 ,6e-04 

8.2e-05 

1 .2e-09 


Figure  7.  Backward  facing  step  in  3D  at  Re  =  1000. 
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Figure  8.  Lift  coefficient  versus  time  for  0  <  Re(t)  <  100. 
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velocity 
1 .4e+00 

1 ,0e+00 

6,8e-01 

3.4e-01 

8.0e-06 


Figure  9.  Flow  around  a  square  cylinder  in  3D  at  Re  =  70  . 
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5  Conclusions 


We  described  a  set  of  numerical  methods  for  approximating  widely  used 
models  of  incompressible  flow.  The  methods  apply  to  complex  geome¬ 
tries  and  unstructured  meshes,  provide  locally  conservative  velocity  fields, 
and  achieve  higher-order  accuracy  in  space  and  time.  Furthermore  the 
multiscale  variational  method  employed  provides  reasonable  accuracy 
for  high  Re  flows  and  has  shown  promise  as  a  hybrid  LES/DNS  method 
(Hoffman  and  Johnson,  2006;  Bazilevs  et  al.,  2007).  The  methods  and 
implementation  were  verified  on  a  set  of  two-  and  three-dimensional 
benchmark  problems.  Several  issues  for  future  work  were  identified: 


•  An  alternative  to  the  quasi-static  subgrid  scales  assumption  in  the 
subgrid  error  approximation  should  be  implemented  for  small  times 
steps. 

•  Futher  work  on  error  estimation  and  startup  heuristics  are  needed 
since  spatial  and  temporal  error  are  tightly  coupled  and  error 

•  Work  vs.  error  studies  should  be  conducted  to  verify  that  the  second- 
order  methods  are  superior  to  first-order  methods. 
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